Phase diagram of imbalanced fermions in optical lattices 
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The zero-temperature phase diagrams of imbalanced fermions in 3D optical lattices are investi- 
gated to evaluate the validity of the Fermi-Hubbard model. It is found that depending on the filling 
factor, s-wave scattering strength and lattice potential, the system may fall into the normal(iV) 
phase, magnetized superfluid(SFju) or phase separation of N and BCS state. By tuning these pa- 
rameters, the superfluidity could be favorable by enhanced effective couplings or suppressed by the 
increased band gap. The phase profiles in the presence of a harmonic trap are also investigated 
under LDA, which show some exotic shell structures compared to those without the optical lattice. 
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In the past few years, great experimental progress has 
been achieved in studying ultracold Fermi gases with 
polarizatio n i 2 i 3 i i ■ With two unequal mixtures of cold 
6 Li atoms in a harmonic trapii 2 -, a clear evidence of phase 
separation with an unpolarized superfluid (BCS) core 
and a normal (N) shell around that has been observed 
in experiment. Theoretically^ ' 8 ' 9 ' 10 ' 11 ' 12 ' 13 , many other 
ground state candidates have been proposed in such 
systems, including magnetized superfluid (SFjf), phase 
separation (PS), and Fulde-Ferrell-Larkin-Ovchinnikov 
(FFLO) state with finite momentum pairing by tun- 
ing the interaction parameter l/kpa s , the polarization 
P = Sn/n or Zeeman field h. Since the optical lattice 
height Vq is also tunable, it is very interesting to study 
its effect on the new phase diagram. For equal mix- 
tures, a second-order quantum phase transition between 
superfluid (SF) and insulating (IN) phases has been ad- 
dressed both experimentally at a critical lattice height V c 
at resonance^ and theoretically 15 '— based on the second- 
order perturbation theory. Besides, there are also works 
on imbalanced fermions in optical lattices focusing on 
INi£, FFLOi£ and SF^i 9 . phases, based on an effective 
Fermi-Hubbard model. 

In this work, starting from the exact lattice spectrum, 
we study the ground state phase diagram of imbalanced 
two species Fermi gases trapped in 3D optical lattices, in 
terms of the total filling factor n, polarization P, s-wave 
scattering length a s and lattice potential Vq. Limited 
by the numerical attainment, the FFLO-type pairing is 
not considered. The total pairing reciprocal lattice mo- 
mentums involved in our simulation are up to the six 
smallest non-zero ones, which turn out to be more and 
more important as Vq increases. Sufficient multiple bands 
have been taken into account to ensure the accuracy es- 
pecially in the strong coupling regime. We demonstrate 
that there are two contradictory effects of Vq on the SF 
phase, depending on the average filling factor n. One is 
the enhanced density of states (DOS) inside each band 
which effectively increases the coupling strength and thus 
is favorable to SF; the other is the broadened band gap 
or discontinuity of DOS which is against SF. One key 
point is that besides tuning a s through the Feshbach res- 
onance(FR), Vq can also be tuned and drive the system 



from weak to strong coupling regime, provided that the 
filling factor is properly fixed. An obvious evidence is the 
emergence of SFm phase for deep optical lattices at par- 
ticular filling regimes, even in far BCS side of FR. We also 
propose that the critical polarization versus total filling 
factor diagram obtained can be used to evaluate the va- 
lidity of the usual Fermi-Hubbard model. The phase pro- 
file in the presence of an external harmonic trap, which is 
more relevant to the practical experiment will be studied 
with local density approximation (LDA) finally. Some 
exotic structures appear, reflecting the uniqueness of the 
optical lattices. 

In a recent experiment**, two hyperfine states of ul- 
tracold 6 Li atoms, \F = 1/2, m F = 1/2) (| |}) and 
\F = l/2,m F = -1/2) (| I)), had been successfully 
loaded to an optical lattice. The low-energy interactions 
are characterized by a single s-wave scattering length a s , 
which can be tuned by FR. Such a system can be well 
described by the one-channel Hamiltonian 

H = J dvY, ^(r)ffo(r)Vv(r) + 



<T=U 

g I dr^(r)^{(r)^(r)^(r), (1) 



where H = £<=*,„,, -^d? /2M + V sin 2 (the; /a); a = 
A/2 is the period of the lattice generated in each direction 
by two oppositely propagating lasers with wavelength A; 
Vq is the lattice height which is usually measured by the 

recoil energy En = 2 Ma 2 ' 9 ls ^ ne renormalized contact 
interaction constant between two species by eliminating 
the unphysical divergence due to the high-momentum 



contribution for fermi gases, - = 



IV J_ 

V Z^q 2e„ ■ 



In the framework of mean-field approach, we expand 
first each field operator in terms of eigenwave functions of 
Ho, ipo-{r) = X)nk^nk(r)V'nk ( T- The Bloch wave functions 
</Wr) = ^Eg a nk (G)e 4 ( k+G )- r and energies e nk are 
obtained from the Schrodinger equation 
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Onk(G') = e„ k a„k(G), (2) 
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where n = {n x ,n y ,n z } = s,p,... indicate the band 
indices; k lie in the first Brillouin zone (BZ) and 
G = 2n /a(l x ,l y ,l z ) is the reciprocal lattice vector. 
The solutions satisfy ^ G a* k (G)a n 'k(G) = S nn > and 
fln,-k(— G) = a* k (G). The standard mean-field treat- 
ment gives 
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Ec^-kf-GMG + Q). Since 

= <5mn and if m ^ n are quite small, in 

the following text we only consider pairing within each 
single band, which means A mn k ~ A n k<5 m n- In such a 
case, the Hamiltonian can be easily diagonalized, and the 
thermodynamic potential is calculated at T = as 



n 
v 



{0( — i? n k+)£'nk+ + ©(— -Enk-)£'nk- 



£nk - M - V ( £ nk - /i) 2 + A^} - 



|A Q | 2 
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with £ nk ± = a/ ( £ nk - m) 2 + A^ T h where /i = (/if + 
/xj.)/2 and h = — /Uj.)/2. From dfi/SAq = and 
iV CT = —d£l/dfi a , we get the gap and density equations 
as 



Aq = j_ v 
9 v ^ 



V 



B„k±>0 



2V( £ nk-M) 2 + A^ k ' 



(6) 
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£ nk± >o V^nk - M) 2 + A^ k ' 

i< E i- E D- (^) 
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Here is the total number of lattice sites, n = (iVj + 
N L )/N L and <5n = (iVf - N L )/N L are the total filling 
factor and the difference, respectively. Hereafter we scale 
all the energies in unit of -Er and the momenta of 2ir/a. 

To make the numerical simulations attainable 
but still retain the essence of the problem, be- 
sides Q = we consider other six non-zero Q: 
(±1, 0, 0), (0, ±1, 0), (0, 0, ±1). Due to = Af± ± Q k and 
the isotropy of 3D cubic lattices, all six non-zero Q share 
the same pairing amplitude Aj. Therefore we get two 
coupled gap equations in terms of A and Aj. For a real- 
istic numerical simulation, we apply a cutoff momentum 



1 2A. | = § (1, 1, 1) in the renormalization and correspond- 
ingly consider the lowest three bands(s,p,c0 in each di- 
rection of lattice spectrum. This truncation allows totally 
n = 54 atoms per site at most, which is well above the 
experimental interest as well as ours in this paper. 




FIG. 1: (color online) Ao, Ai and Ai/Ao(see inset) vs lattice 
potentials Vo for equal mixtures. a/a s — —3. The averaged 
filling is fixed to be n = 1. 

Before turning to the phase diagram of imbalanced 
system, first we analyze the necessity of involving non- 
zero Q in gap equations for equal mixtures. Fig. [T] 
shows Ao,Ai and their ratio as a function of lattice po- 
tential Vq at half filling n = 1. It is shown that the 
Q ^ pairing becomes more and more important as 
Vo increases. This effect can be understood as follows. 
Taking a very shallow ID lattice for example, non-zero 
m51 with \Q\ = 2,4,6... and 1,3,5... only exist around 
kinetic energy-degenerate points k = and k = ±1/2 
respectively, which contribute little to gap-equations and 
therefore produce a negligible Ai. In the limit of Vq = 0, 
these non-zero M® k exactly cancel with each other in 
gap equations and finally only Q = pairing survives. 
However as Vq increases, the eigen- vector a n k{G) evolves 
such that the area of non-zero MR expands from three 
discrete points in first BZ to considerable regions around 
them, leading to an increasing Aq with Q ^ 0. Since our 
interest is still within s-band, the \Q\ = 1 pairings take 
a leading role among all non-zero ones, which is verified 
both numerically and analytically from a perturbation 
theory for shallow lattices. This is why we just take into 
account six smallest non-zero Q in 3D case for not-so- 
deep lattices. The consideration of non-zero Q-pairing 
would produce a much stronger superfluidity especially 
for deep lattices, which can also be seen from the com- 
parison of the previous two works^. 

The ground state phase diagram in Fig. [5] is de- 
termined as follows. We first compare flBCsiPih = 
E min ,A = {A ,Ai}) with with Q N (fi,h = E min ,A = 
{0}), with Aq/j obtained for unpolarized BCS state and 
E mm = Min„ fe ( v / (e nfc - /i) 2 + Af lk ) its lowest excitation 
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energy. If flscs > then the first-order phase transi- 
tion point h c {< E min ) is given by 

Sl B cs(»,h c ,A = {A ,A 1 }) = Sl N (n,h c ,A = {Q})[8) 
n N (p,,h c ) = n. (9) 

P c = 5riN(fi, h c )/n represents a critical point when PS 
is entirely composed by N phase. Note that in this case 
the polarized SF or Sarma phased is unstable due to the 
negative superfluid density- 11 . If Qbcs < then there 
should be a stable SFm interpolating between BCS and 
N phase. In free space at the SFm-N second-order transi- 
tion point, N denotes a fully polarized normal state with 
P c = l 13 . Correspondingly in optical lattices, we obtain 
P c = 1 at n <C 1 and |n — 2|/n at n ~ 2, as shown by red 
solid circles in Fig. 
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FIG. 2: (color online) Zero temperature phase diagram as a 
function of polarization P c — Sn/n, total filling factor n and 
lattice height Vo- a/a s = —3. The red dashed(blue solid) 
lines show P c evolves with Vo(n) for fixed n(Vo). All the lines 
above denote the PS — N boundaries, expect for red solid 
circles separating SFm and N phase instead. 

We analyze that P c — Vo curves reveal two contradic- 
tory effects of increasing Vo to SF depending on the filling 
factor n. As shown in Fig. [31 for n < 1 increasing Vo 
will flatten each band and enhance DOS (almost inversely 
proportional to the band width); while at n ~ 2, increas- 
ing Vo produces an entirely opposite effect due to the 
enlarged band gap. According to the standard BCS the- 
ory, the DOS at the Fermi surface dramatically affects 
the strength of SF, as is also reflected by such contra- 
dictory effects. When Vo = 3Er, P c increases to unity 
at small n but reduces to zero at n = 2, denoting the 
IN phase with n-f = n± = 1. For n e (1,2), P c initially 
drops down and then goes up, indicating the competition 
between the above these two effects. Here the lattice en- 
hancement of P c at n < 1 is similar to the enhancement 
of T c for equal mixtures in weak coupling limi t 21 ' 22 . 

Next we turn to P c — n diagram for fixed Vo. As 
is well known in free space, a first-order BCS to N 
phase transition takes place in weak coupling limit at 
h c = A= and P c — with the gap amplitude 
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FIG. 3: Density of state(DOS) at the Fermi surface versus 
filling factor n for single-spin atoms in 3D free space and in 
lattices with Vq/Er = l(no band gap) and Vq/Er = 3(with 
band gap). Inset is DOS for non-interacting Hubbard model. 
The dotted lines therein denote two peaks of DOS at (/i = 
t, n p = 0.213) and (3t — fj.,1 — n p ). 

Aq = §zEpexp(— 2fcj *j a | ) and the interaction parame- 
ter ri = , 1 = — (37r 2 n)~3. p will increase with n all 

' kpa s a s v > L 

along from weak coupling limit(^ — > -co) to the unitary 
limits — > 0). Within an optical lattice, however, the 
P c — n curve would be dramatically modified. In weak 
coupling limit with small Vo, the curve basically follows 
as that of DOS in Fig. [3l with a dip at rid ~ 2 and 
correspondingly a peak at n p . As Vo increases, n v grad- 
ually moves to the left and finally vanishes at n p = 0, 
and finally SFm state emerges at n <C 1 or n ~ 2. Dif- 
ferent from the SFm studied by DMFT method under 
tight-binding modet^, the phase shown here is purely 
due to the enhanced effective coupling by lattices. In 
this limit, two fermions are likely to form a molecule, 
and the BCS equation directly reduces to a Schrodinger 
equation for a single bound pai r 21 ' 23 . It is expected that 
as Vo increases, the SF M phase would extend to a larger 
or even the whole density region. Actually, the physics 
at n <C 1 and n ~ 2 can be related to each other via 
particle-hole symmetry. The symmetry is essentially ob- 
vious within the background of Fermi-Hubbard model, 
£k = 0.5i — cos(fcid)). Since {n,^) and (2 — n, 3< — ii) 
share the same {A, h, Sn, fi} and thus the same critical 
Sn c , the critical polarization P c (n) for n < 8 follows as 

f 1, n<l 
P c ={ JnrHL, 1< n < 5 (10) 
I 5<n<8 

We also compute the phase diagram at other s-wave 
couplings with fixed Vb as shown in Fig. 2J Different 
from Fig. [2 it shows that the increasing a/a s always en- 
hance SF and improve P c regardless of filling factors. At 
sufficiently strong coupling close to unitary, the particle- 
hole symmetry in each band breaks down since it is en- 
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ergetically favorable for particles in s-band to overcome 
the band-gap and form cooper pairs even at n = 2. In 
this case the multi-band effect should be taken into ac- 
count. This is why SFjvf only turns up at n <C 1 but not 
at n ~ 2 in unitary limit, shown as blue circles in Fig. [4] 
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FIG. 4: (color online) P c versus n diagram at different cou- 
plings. Vo = 3Er. All lines denote PS — N boundaries, except 
that the green and blue circles show SFm — N boundaries. 
The black dashed-dot line marks the upper limit of peak po- 
sition based on Fermi- Hubbard model(see text). 

We emphasize that the P c — n diagram in weak cou- 
pling limit can be used to evaluate the validity of tight- 
binding approximation(TBA) usually applied to the lat- 
tices. For Hubbard model, the DOS shows two peaks 
symmetric around half filling(see inset of Fig. [3]), due 
to the van Hove singularity at (fj, = t,n = 0.213) and 
(/j = 2t,n = 0.787). We also verified numerically that 
the peak position of P c at different couplings is never 
greater than 0.426 for arbitrary interactions \U\/t, twice 
as that for the first peak in DOS. This universal prop- 



erty could be used to justify the validity of TBA to re- 
alistic lattices. Apparently from Fig. SI we see that the 
TBA is not applicable to Vq = 3Er, since at really weak 
coupling(a/a s = —9) the peak position n p ~ 0.6 > 0.426. 
The disagreement here indicates the deviation of the two 
lattice spectrum, and thus the necessity of adopting exact 
lattice spectrum for not-so-deep lattices. 

Finally, it is also useful to consider the phase profile rel- 
evant to realistic experiments with an external harmonic 
potential V(r). Under LDA, the system is assumed to be 
locally homogeneous with an averaged chemical potential 
/i(r) = (^ot + A*oi)/2 — V(r) and position-independent 
difference h = (/i j — /xoj.)/2. The phase at position r 
is determined by the local (fi(r),h), which is also self- 
consistently related to the total particle numbers, s-wave 
interaction and the lattice potential. Here we give sev- 
eral typical phase profiles with the filling factor in trap 
center larger than 2. For relatively shallow lattices and 
very weak s-wave interactions, a typical one from the 
trap center to the edge is: BCS-PN-IN-PN-BCS-PN- 
FN (PN/FN: partially/fully polarized Normal). Start- 
ing from this profile, if Vq increases, PN is very likely 
to be replaced by SFm at positions where n(r) ~ 2 or 
n(r) -C 1, while IN still survive in a certain region; But if 
s-wave interaction increases, then IN shrinks gradually, 
two BCS regimes merge together and PN gives rise to 
SF M , with only three phases left finally: BCS-SF M -FN. 
In the latter case, much higher bands with continuous 
spectrum would be occupied, which makes the situation 
very similar to free space and the lattice effect is not 
obvious. 
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